eea2020_tp2_dario_hruszecki

# Carga de librerías
library(tidyverse)

## ── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ dplyr   1.0.0
## ✓ tidyr   1.1.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0

## ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(tidymodels)

## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidymodels 0.1.1 ──

## ✓ broom     0.7.0      ✓ recipes   0.1.13
## ✓ dials     0.0.9      ✓ rsample   0.0.8 
## ✓ infer     0.5.3      ✓ tune      0.1.1 
## ✓ modeldata 0.1.0      ✓ workflows 0.2.1 
## ✓ parsnip   0.1.4      ✓ yardstick 0.0.7

## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidymodels_conflicts() ──
## x scales::discard() masks purrr::discard()
## x dplyr::filter()   masks stats::filter()
## x recipes::fixed()  masks stringr::fixed()
## x dplyr::lag()      masks stats::lag()
## x yardstick::spec() masks readr::spec()
## x recipes::step()   masks stats::step()

library(GGally)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

library(corrr)
library(knitr)
library(ggrepel)
library(kableExtra)

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

library(here)

## here() starts at /home/daro/cdatos/TPs/EEA/TP2

EEA 2020 - TP2 - Darío Hruszecki

DATOS

El dataset de precios de inmuebles proviene de Properati. El mismo ya fue filtrado por los docentes y a su vez se encuentra particionado en subconjuntos de training y test.

Carga partición de entrenamiento

train_ds <- read_csv(here("/ds/ar_properties_train.csv"))

## Parsed with column specification:
## cols(
##   id = col_character(),
##   l3 = col_character(),
##   rooms = col_double(),
##   bathrooms = col_double(),
##   surface_total = col_double(),
##   surface_covered = col_double(),
##   price = col_double(),
##   property_type = col_character()
## )

Describimos el dataset

Previsualizamos los datos

head(train_ds)

## # A tibble: 6 x 8
##   id    l3    rooms bathrooms surface_total surface_covered  price property_type
##   <chr> <chr> <dbl>     <dbl>         <dbl>           <dbl>  <dbl> <chr>        
## 1 r22O… Alma…     1         1            40              37  92294 Departamento 
## 2 atZQ… Alma…     1         1            49              44 115000 Departamento 
## 3 Tur/… Alma…     1         1            40              37  88900 Departamento 
## 4 PNtN… Alma…     1         1            40              37  88798 Departamento 
## 5 Lhg4… Alma…     1         1            49              44 110975 Departamento 
## 6 oeLh… Alma…     1         1            40              37  99000 Departamento

Vemos las dimensiones del DS y como se distribuyen los valores de las variables

dim_desc(train_ds)

## [1] "[32,132 x 8]"

Vemos como se distribuyen los valores de las variables

summary(train_ds)

##       id                 l3                rooms         bathrooms    
##  Length:32132       Length:32132       Min.   :1.000   Min.   :1.000  
##  Class :character   Class :character   1st Qu.:2.000   1st Qu.:1.000  
##  Mode  :character   Mode  :character   Median :3.000   Median :1.000  
##                                        Mean   :2.722   Mean   :1.427  
##                                        3rd Qu.:3.000   3rd Qu.:2.000  
##                                        Max.   :8.000   Max.   :5.000  
##  surface_total surface_covered      price        property_type     
##  Min.   : 28   Min.   : 26.00   Min.   : 69500   Length:32132      
##  1st Qu.: 46   1st Qu.: 42.00   1st Qu.:120000   Class :character  
##  Median : 65   Median : 58.00   Median :169000   Mode  :character  
##  Mean   : 80   Mean   : 70.61   Mean   :214991                     
##  3rd Qu.: 99   3rd Qu.: 85.00   3rd Qu.:260000                     
##  Max.   :320   Max.   :265.00   Max.   :950000

Valores únicos para la variable alfanumérica l3

unique(train_ds$l3)

##  [1] "Almagro"              "Palermo"              "Caballito"           
##  [4] "Villa Crespo"         "Villa Urquiza"        "Barracas"            
##  [7] "Nuñez"                "Belgrano"             "Floresta"            
## [10] "Recoleta"             "Colegiales"           "Parque Chas"         
## [13] "Barrio Norte"         "Villa Devoto"         "Villa Ortuzar"       
## [16] "Velez Sarsfield"      "Parque Chacabuco"     "Villa Pueyrredón"    
## [19] "Villa Real"           "Once"                 "Flores"              
## [22] "San Telmo"            "Las Cañitas"          "Villa Santa Rita"    
## [25] "Villa del Parque"     "Retiro"               "Parque Centenario"   
## [28] "Chacarita"            "Abasto"               "San Cristobal"       
## [31] "Liniers"              "Balvanera"            "Boedo"               
## [34] "Villa General Mitre"  "Saavedra"             "Parque Patricios"    
## [37] "Boca"                 "Paternal"             "Monserrat"           
## [40] "San Nicolás"          "Parque Avellaneda"    "Centro / Microcentro"
## [43] "Puerto Madero"        "Congreso"             "Monte Castro"        
## [46] "Mataderos"            "Coghlan"              "Versalles"           
## [49] "Villa Lugano"         "Catalinas"            "Villa Luro"          
## [52] "Constitución"         "Pompeya"              "Tribunales"          
## [55] "Agronomía"            "Villa Soldati"        "Villa Riachuelo"

Valores únicos para la variable alfanumérica property_type

unique(train_ds$property_type)

## [1] "Departamento" "PH"           "Casa"

Verificamos si hay valores NA en las variables

apply(train_ds, 2, function(x) any(is.na(x)))

##              id              l3           rooms       bathrooms   surface_total 
##           FALSE           FALSE           FALSE           FALSE           FALSE 
## surface_covered           price   property_type 
##           FALSE           FALSE           FALSE

table(train_ds$property_type)

## 
##         Casa Departamento           PH 
##          809        28203         3120

tail(table(train_ds$l3)[], 10)

## 
## Villa General Mitre        Villa Lugano          Villa Luro       Villa Ortuzar 
##                  99                 113                 264                 109 
##    Villa Pueyrredón          Villa Real     Villa Riachuelo    Villa Santa Rita 
##                 285                  55                  13                 128 
##       Villa Soldati       Villa Urquiza 
##                   9                1319

Observaciones sobre los datos

• El dataset contiene 8 variables y 32132 registros
• De las 8 variables, 5 son núnericas y 3 son alfanuméricas
• De las variables alfanuméricas tenemos el barrio (l3) y el tipo de propiedad (property_type). La tercera es un UUID para la propiedad (id)
• En base a los valores únicos encontrados para l3 podemos inferir que las propiedades pertenecen a CABA en su totalidad.
• De las númericas, además del precio (price) tenemos variables que describen la propiedad: cantidad de habitaciones (rooms), cantidad de baños (bathrooms), supervicie total (surface_total) y superficie cubierta (surface_covered)
• No se observan valores NA en ninguna de las columnas

CONSIGNA

Modelo de Regresion Lineal Multiple

Generación del Modelo con todas las covariables

Vamos a generar un modelo que utilice todas las covariables presentes en nuestro dataset de training para buscar predecir el precio.

modelo_todas_covariables <- lm(
  price ~ rooms + bathrooms + surface_total + surface_covered + property_type + l3, 
  data = train_ds
  )
tidy_mtc <- tidy(modelo_todas_covariables, conf.int = TRUE)
tidy_mtc

## # A tibble: 63 x 7
##    term                estimate std.error statistic   p.value conf.low conf.high
##    <chr>                  <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
##  1 (Intercept)         -109144.    5784.    -18.9   5.28e- 79 -120480.   -97808.
##  2 rooms                 -4360.     530.     -8.23  2.01e- 16   -5399.    -3321.
##  3 bathrooms             34367.     768.     44.7   0.          32862.    35873.
##  4 surface_total           891.      28.2    31.5   3.84e-215     835.      946.
##  5 surface_covered        1498.      34.3    43.6   0.           1431.     1565.
##  6 property_typeDepar…   91486.    2660.     34.4   1.35e-254   86272.    96700.
##  7 property_typePH       46220.    2752.     16.8   5.06e- 63   40826.    51615.
##  8 l3Agronomía            1241.   10652.      0.117 9.07e-  1  -19638.    22120.
##  9 l3Almagro             -3480.    5182.     -0.672 5.02e-  1  -13636.     6676.
## 10 l3Balvanera          -25835.    5491.     -4.70  2.55e-  6  -36598.   -15072.
## # … with 53 more rows

Análisis del Modelo con todas las covariables

Significado de los coeficientes estimados

• El valor de la ordenada al origen (-109.143,78) es el valor del precio esperado para una propiedad sin ambientes o superficie y que no pertenece a ningún tipo (Casa, PH o Departamento) y que no se encuentra en ninguna zona (es una inferencia más teórica de nuestro modelo dado que no es representativo en la realidad).

• El coeficiente estimado de ambientes es -4359,73. Esto implica que manteniendo constantes todas nuestras otras covariables, un aumento de un ambiente a nuestra propiedad corresponde a una quita de $4359, en promedio en el precio de una propiedad.

• El coeficiente estimado de baños es 34.367,21. Esto implica que manteniendo constantes todas nuestras otras covariables, un aumento de un baño a nuestra propiedad corresponde a una aumento de $34.367, en promedio en el precio de una propiedad.

• El coeficiente estimado de superficie total es 890,56. Esto implica que manteniendo constantes todas nuestras otras covariables, un aumento en una unidad de superficie (m2) a nuestra propiedad corresponde a una aumento de $890, en promedio en el precio de una propiedad.

• El coeficiente estimado de superficie cubierta es 1.497,93. Esto implica que manteniendo constantes todas nuestras otras covariables, un aumento en una unidad de superficie cubierta (m2) a nuestra propiedad corresponde a una aumento de $1.498, en promedio en el precio de una propiedad.

• Los coeficientes de tipo de propiedad (Departamento y PH) indican cuanto aumenta la función de respuesta (precio) para un departamento o un ph en comparación a una casa (categoría basal). Sus valores de β respectivos son para los departamentos 91.485,94 y para los PHs 46.220,04

• Finalmente tenemos 56 coeficientes para la variable l3 que representan en cuanto aumenta o se reduce el precio de una propiedad dependiendo de la zona específica a la cual pertenece.

Significatividad de variables dummies

Observamos en nuestro modelo dos tipos distintos de variables dummies. Aquellas que responden a la variable property_type y aquellas que responden a la variable l3.

Podemos observar que los p-valores asociados a los coeficientes de la variable property_type, resultan estadísticamente significativos. En cuanto a aquellos p-valores asociados a los coeficientes de la variable l3 nos encontramos con que en algunos casos resultan estadísticamente significativos (l3Balvanera con un p-valor de 0,000002) pero en otros casos no lo es (l3Agronomía con un p-valor de 0,9).

Además podemos observar que los intervalos de confianza para las variables dummy que responden a property_type no contienen al 0 y en lineas generales, aquellos que responden a la variable l3 si lo contienen.

Vamos a aplicar un Test F (test de ignificatividad global) para medir la significatividad conjunta de las variables categóricas en nuestro modelo.

tidy(
  anova(modelo_todas_covariables)
)

## # A tibble: 7 x 6
##   term               df   sumsq  meansq statistic p.value
##   <chr>           <int>   <dbl>   <dbl>     <dbl>   <dbl>
## 1 rooms               1 2.18e14 2.18e14    48554.       0
## 2 bathrooms           1 1.09e14 1.09e14    24391.       0
## 3 surface_total       1 7.74e13 7.74e13    17276.       0
## 4 surface_covered     1 1.73e13 1.73e13     3863.       0
## 5 property_type       2 2.03e13 1.01e13     2265.       0
## 6 l3                 56 5.90e13 1.05e12      235.       0
## 7 Residuals       32069 1.44e14 4.48e 9       NA       NA

La tabla de ANOVA muestra que tanto la variable categorica property_type como l3 en su conjunto resultan estadísticamente significativas (p-valor < 0.05). Por ende buscaremos determinar cuales de las comparaciones entre grupos son estadísticamente significativas.

Evaluación del Modelo

Graficamos los coeficientes estimados.

ggplot(
  tidy_mtc, 
  aes(
    estimate, 
    term, 
    xmin = conf.low, 
    xmax = conf.high, 
    height = 0
  )
) +
  geom_point(
    color = "forestgreen", 
    size=2
  ) +
  geom_vline(
    xintercept = 0, 
    lty = 4, 
    color = "black"
  ) +
  geom_errorbarh(
    color = "forestgreen", 
    size=1
  ) + 
  theme_bw() +
  labs(
    y = "Coeficientes β", 
    x = "Estimación"
  )

De analizar el grafico de los coeficientes podemos visualizar los intervalos de confianza de las distintas variables. Claramente podemos observar como hay muchas variables dummy asociadas a la variable l3 que incluyen el 0 dentro de sus intervalos de confianza y por ende podemos asumir que estas variables no son significativas para explicar el precio de las propiedades.

Y ahora analizamos las medidas del resumen global del modelo.

glance(modelo_todas_covariables)

## # A tibble: 1 x 12
##   r.squared adj.r.squared  sigma statistic p.value    df  logLik    AIC    BIC
##       <dbl>         <dbl>  <dbl>     <dbl>   <dbl> <dbl>   <dbl>  <dbl>  <dbl>
## 1     0.777         0.777 66933.     1803.       0    62 -4.03e5 8.05e5 8.06e5
## # … with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

Analizando el valor del \(Rˆ2\) ajustado del modelo, podemos ver que el mismo parece explicar un buen porcentaje de variabilidad del conjunto.

### Generación de Modelo sin l3

Ahora vamos a generar un modelo sin la variable l3 para comparar y analizar el impacto de esta variable en nuestro modelo previo.

modelo_sin_l3 <- lm(
  price ~ rooms + bathrooms + surface_total + surface_covered + property_type, 
  data = train_ds
)
tidy_msl3 <- tidy(modelo_sin_l3, conf.int = TRUE)
tidy_msl3

## # A tibble: 7 x 7
##   term                 estimate std.error statistic   p.value conf.low conf.high
##   <chr>                   <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)          -128676.    3333.      -38.6 1.07e-318 -135210.  -122143.
## 2 rooms                 -13657.     617.      -22.1 8.68e-108  -14867.   -12448.
## 3 bathrooms              42564.     899.       47.3 0.          40801.    44327.
## 4 surface_total            856.      33.0      25.9 1.40e-146     791.      920.
## 5 surface_covered         1819.      40.0      45.5 0.           1740.     1897.
## 6 property_typeDepart…  133033.    3049.       43.6 0.         127056.   139010.
## 7 property_typePH        66457.    3233.       20.6 2.59e- 93   60121.    72793.

Podemos observar como en nuestro nuevo modelo, los p-valores de todos los coeficientes asociados a nuestras variables resultan estadísticamente significativos (p-valor < 0.05). Además ninguno de los intervalos de confianza de nuestros coeficientes incluyen el 0.

Vamos a aplicar un Test F (test de ignificatividad global) para medir la significatividad conjunta de las variables categóricas en nuestro modelo.

tidy(
  anova(modelo_sin_l3)
)

## # A tibble: 6 x 6
##   term               df   sumsq  meansq statistic p.value
##   <chr>           <int>   <dbl>   <dbl>     <dbl>   <dbl>
## 1 rooms               1 2.18e14 2.18e14    34486.       0
## 2 bathrooms           1 1.09e14 1.09e14    17324.       0
## 3 surface_total       1 7.74e13 7.74e13    12271.       0
## 4 surface_covered     1 1.73e13 1.73e13     2744.       0
## 5 property_type       2 2.03e13 1.01e13     1609.       0
## 6 Residuals       32125 2.03e14 6.31e 9       NA       NA

Los resultados del ANOVA muestran que todas las variables utilizadas son estadisticamente significativas (respaldando lo visto en el anterior modelo).

Y ahora analizamos comparativamente las medidas del resumen global de nuestros 2 modelos.

modelos = list(
  modelo_todas_covariables = modelo_todas_covariables,
  modelo_sin_l3 = modelo_sin_l3
)

purrr::map_df(modelos, broom::glance, .id = "model")

## # A tibble: 2 x 13
##   model r.squared adj.r.squared  sigma statistic p.value    df  logLik    AIC
##   <chr>     <dbl>         <dbl>  <dbl>     <dbl>   <dbl> <dbl>   <dbl>  <dbl>
## 1 mode…     0.777         0.777 66933.     1803.       0    62 -4.03e5 8.05e5
## 2 mode…     0.686         0.686 79419.    11674.       0     6 -4.08e5 8.16e5
## # … with 4 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>,
## #   nobs <int>

Como se puede observar, el modelo_todas_covariables es el modelo que mejor explica la variabilidad en nuestro conjunto (evaluamos el coeficiente \(Rˆ2\) ajustado para ambos).

Creación de Variables

Armado de nuevo dataset

En nuestro analisis previo, observamos como algunos barrios (variable l3) eran significativos y otros no. Por ende, vamos a proponer la generación de una nueva variable barrios que nos permite agrupar mas eficientemente nuestros datos.

Para hacerla, vamos a generar otra variable precio_m2 que contiene el valor de price sobre el total de metros cuadrados de superficie (surface_total). Luego, sobre esta variable, agruparemos los datos en tres categorías (variable barrios): - precio_bajo: aquellas propiedades con valores iguales o inferiores al Q1 de nuestra variable precio_m2 - precio_medio: aquellas propiedades con valores superiores al Q1 y con valores inferiores o iguales al Q3 de nuestra variable precio_m2 - precio_alto: aquellas propiedades con valores superiores al Q3 de nuestra variable precio_m2

train_ds = train_ds %>%
  mutate(
    precio_m2 = price / surface_total,
    barrios = factor(
      case_when(
        precio_m2 <= quantile(precio_m2)[2] ~ "precio_bajo",
        precio_m2 > quantile(precio_m2)[2] & precio_m2 <= quantile(precio_m2)[4] ~ "precio_medio",
        precio_m2 > quantile(precio_m2)[4] ~ "precio_alto"
      )
    )
  )

train_ds

## # A tibble: 32,132 x 10
##    id    l3    rooms bathrooms surface_total surface_covered  price
##    <chr> <chr> <dbl>     <dbl>         <dbl>           <dbl>  <dbl>
##  1 r22O… Alma…     1         1            40              37  92294
##  2 atZQ… Alma…     1         1            49              44 115000
##  3 Tur/… Alma…     1         1            40              37  88900
##  4 PNtN… Alma…     1         1            40              37  88798
##  5 Lhg4… Alma…     1         1            49              44 110975
##  6 oeLh… Alma…     1         1            40              37  99000
##  7 Fmtr… Alma…     1         1            40              37  96984
##  8 QohS… Pale…     1         1            40              34  99500
##  9 WEJo… Pale…     1         1            36              33 121600
## 10 ByY9… Pale…     3         2            90              90 285000
## # … with 32,122 more rows, and 3 more variables: property_type <chr>,
## #   precio_m2 <dbl>, barrios <fct>

Generación de un nuevo modelo

Ahora vamos a generar un nuevo modelo incluyendo nuestras nueva covariables barrios.

modelo_barrios <- lm(
  price ~ rooms + bathrooms + surface_total + surface_covered + property_type + barrios , 
  data = train_ds
) 
tidy_mb <- tidy(modelo_barrios, conf.int = TRUE)
tidy_mb

## # A tibble: 9 x 7
##   term                 estimate std.error statistic   p.value conf.low conf.high
##   <chr>                   <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)            -3082.    2574.      -1.20 2.31e-  1   -8127.     1962.
## 2 rooms                  -3831.     458.      -8.37 6.08e- 17   -4728.    -2934.
## 3 bathrooms              21781.     674.      32.3  2.58e-225   20459.    23102.
## 4 surface_total           1561.      24.9     62.8  0.           1512.     1610.
## 5 surface_covered         1103.      29.8     37.0  1.22e-293    1044.     1161.
## 6 property_typeDepart…   87739.    2266.      38.7  1.02e-320   83299.    92180.
## 7 property_typePH        59258.    2379.      24.9  1.14e-135   54596.    63921.
## 8 barriosprecio_bajo   -165257.    1017.    -163.   0.        -167250.  -163264.
## 9 barriosprecio_medio   -93607.     810.    -116.   0.         -95195.   -92020.

Podemos apreciar como los p-valores de los coeficientes resultan significativos. Sin embargo el valor de mi coeficiente \(B_0\) parece no ser significativo (p-valor > 0.05). Más aún, podemos comprobar que su intervalo de confianza incluye el 0.

Evaluación Modelo Barrios

Comparamos nuestros 3 modelos.

modelos = list(
  modelo_todas_covariables = modelo_todas_covariables,
  modelo_sin_l3 = modelo_sin_l3,
  modelo_barrios = modelo_barrios
)

purrr::map_df(modelos, broom::glance, .id = "model") %>%
  arrange(adj.r.squared)

## # A tibble: 3 x 13
##   model r.squared adj.r.squared  sigma statistic p.value    df  logLik    AIC
##   <chr>     <dbl>         <dbl>  <dbl>     <dbl>   <dbl> <dbl>   <dbl>  <dbl>
## 1 mode…     0.686         0.686 79419.    11674.       0     6 -4.08e5 8.16e5
## 2 mode…     0.777         0.777 66933.     1803.       0    62 -4.03e5 8.05e5
## 3 mode…     0.830         0.830 58434.    19576.       0     8 -3.98e5 7.97e5
## # … with 4 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>,
## #   nobs <int>

Podemos observar como el modelo que mejor explica la variabilidad, \(Rˆ2\) ajustado, es nuestro modelo_barrios seguido en orden por modelo_todas_covariables y modelo_sin_l3.

Más aún, como el modelo_barrios tiene coeficientes que son estadísticamente significativos, resulta más útil para el análisis del precio de nuestras propiedades (a pesar que la información de los valores de l3 podría ser más util en un análisis mas profundo).

Generación Modelo Supercicie Descubierta

Dada la fuerte correlación entre las variables surface_total y surface_covered, vamos a generar una nueva variable llamada sup_descubierta que representa la diferencia entre la superficie total y la cubierta.

train_ds = train_ds %>%
  mutate(sup_descubierta = surface_total - surface_covered)

train_ds 

## # A tibble: 32,132 x 11
##    id    l3    rooms bathrooms surface_total surface_covered  price
##    <chr> <chr> <dbl>     <dbl>         <dbl>           <dbl>  <dbl>
##  1 r22O… Alma…     1         1            40              37  92294
##  2 atZQ… Alma…     1         1            49              44 115000
##  3 Tur/… Alma…     1         1            40              37  88900
##  4 PNtN… Alma…     1         1            40              37  88798
##  5 Lhg4… Alma…     1         1            49              44 110975
##  6 oeLh… Alma…     1         1            40              37  99000
##  7 Fmtr… Alma…     1         1            40              37  96984
##  8 QohS… Pale…     1         1            40              34  99500
##  9 WEJo… Pale…     1         1            36              33 121600
## 10 ByY9… Pale…     3         2            90              90 285000
## # … with 32,122 more rows, and 4 more variables: property_type <chr>,
## #   precio_m2 <dbl>, barrios <fct>, sup_descubierta <dbl>

Ahora con nuestra nueva variable, vamos a generar un nuevo modelo que incluya la variable sup_descubierta en reemplazo de la variable surface_total.

modelo_sup_descubierta <- lm(
  price ~ rooms + bathrooms + sup_descubierta + surface_covered + property_type + precio_m2 + barrios , 
  data = train_ds
) 
tidy_msd <- tidy(modelo_sup_descubierta, conf.int = TRUE)
tidy_msd

## # A tibble: 10 x 7
##    term                estimate std.error statistic   p.value conf.low conf.high
##    <chr>                  <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
##  1 (Intercept)        -343392.   2485.      -138.   0.         -3.48e5 -338521. 
##  2 rooms                 1064.    313.         3.41 6.61e-  4   4.52e2    1677. 
##  3 bathrooms            12549.    461.        27.2  3.96e-161   1.16e4   13453. 
##  4 sup_descubierta       1970.     17.1      115.   0.          1.94e3    2004. 
##  5 surface_covered       2611.     10.5      249.   0.          2.59e3    2632. 
##  6 property_typeDepa…   68762.   1545.        44.5  0.          6.57e4   71790. 
##  7 property_typePH      51785.   1619.        32.0  6.01e-221   4.86e4   54959. 
##  8 precio_m2               92.2     0.478    193.   0.          9.13e1      93.2
##  9 barriosprecio_bajo   19276.   1180.        16.3  1.00e- 59   1.70e4   21590. 
## 10 barriosprecio_med…   18396.    800.        23.0  5.87e-116   1.68e4   19965.

Podemos observar como los p-valores de todos nuestros coeficientes resultan significativos en nuestro nuevo modelo.

tidy(
  anova(modelo_sup_descubierta)
)

## # A tibble: 8 x 6
##   term               df   sumsq  meansq statistic    p.value
##   <chr>           <int>   <dbl>   <dbl>     <dbl>      <dbl>
## 1 rooms               1 2.18e14 2.18e14   137576.  0.       
## 2 bathrooms           1 1.09e14 1.09e14    69112.  0.       
## 3 sup_descubierta     1 3.81e12 3.81e12     2407.  0.       
## 4 surface_covered     1 9.09e13 9.09e13    57490.  0.       
## 5 property_type       2 2.03e13 1.01e13     6417.  0.       
## 6 precio_m2           1 1.51e14 1.51e14    95488.  0.       
## 7 barrios             2 8.63e11 4.32e11      273.  2.78e-118
## 8 Residuals       32122 5.08e13 1.58e 9       NA  NA

Y al ejecutar el Test F, validamos el mismo.

Ahora comparamos todos los modelos que generamos.

modelos = list(
  modelo_todas_covariables = modelo_todas_covariables,
  modelo_sin_l3 = modelo_sin_l3,
  modelo_barrios = modelo_barrios,
  modelo_sup_descubierta = modelo_sup_descubierta
)

purrr::map_df(modelos, broom::glance, .id = "model") %>%
  arrange(adj.r.squared)

## # A tibble: 4 x 13
##   model r.squared adj.r.squared  sigma statistic p.value    df  logLik    AIC
##   <chr>     <dbl>         <dbl>  <dbl>     <dbl>   <dbl> <dbl>   <dbl>  <dbl>
## 1 mode…     0.686         0.686 79419.    11674.       0     6 -4.08e5 8.16e5
## 2 mode…     0.777         0.777 66933.     1803.       0    62 -4.03e5 8.05e5
## 3 mode…     0.830         0.830 58434.    19576.       0     8 -3.98e5 7.97e5
## 4 mode…     0.921         0.921 39763.    41717.       0     9 -3.86e5 7.72e5
## # … with 4 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>,
## #   nobs <int>

Podemos observar como nuestro nuevo modelo modelo_sup_descubierta es el que mejor explica la variabilidad.

Diagnóstico del Modelo

Vamos a evaluar nuestro modelo lineal con la variable sup_descubierta

plot(modelo_sup_descubierta)

- Residuos vs valores predichos: parece existir una cierta estructura en los datos. La linealidad parece que no se respeta. Vemos la existencía de heterocedasticidad a medida que aumentan las predicciones. Se resaltan 2 observaciones. - Normal QQ plot: los extremos no se ajustan a la distribución teórica, tanto inferior izquierdo como superior derecho. - Residual vs leverage:Existen cuatro puntos con un leverage bastante alto. Y se identifican 2 puntos como posibles outliers. .

• Diagnóstico del modelo: El modelo creado no cumple con los supuestos del modelo lineal. No parece ser concluyente respecto de la homocedasticidad, definitivamente existe una falta de normalidad y hay presencia de observaciones de alto leverage

Modelo Log(price)

Vamos a generar un modelo para predecir el logaritmo natural de la variable price \[ log(price) = \beta_0 + \beta_1log(rooms) + \beta_2log(bathrooms) + \beta_3log(surface\_covered) + \beta_4property\_type + \beta_5barrio + \beta_6surface\_patio \]

Para esto vamos a generar nuestros nuevas variables

modelo_log = lm(
  log(price) ~ log(rooms) + log(bathrooms) + log(surface_covered) + property_type + barrios + sup_descubierta,
  data = train_ds
)
tidy_ml <- tidy(modelo_log, conf.int = TRUE)
tidy_ml

## # A tibble: 9 x 7
##   term                 estimate std.error statistic   p.value conf.low conf.high
##   <chr>                   <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)           8.53    0.0155       549.   0.         8.50      8.56   
## 2 log(rooms)           -0.0211  0.00315       -6.68 2.36e- 11 -0.0273   -0.0149 
## 3 log(bathrooms)        0.0694  0.00324       21.4  2.68e-101  0.0630    0.0757 
## 4 log(surface_covered)  0.901   0.00371      243.   0.         0.894     0.908  
## 5 property_typeDepart…  0.171   0.00611       27.9  1.00e-169  0.159     0.183  
## 6 property_typePH       0.0921  0.00644       14.3  3.61e- 46  0.0794    0.105  
## 7 barriosprecio_bajo   -0.708   0.00278     -254.   0.        -0.713    -0.702  
## 8 barriosprecio_medio  -0.359   0.00222     -162.   0.        -0.364    -0.355  
## 9 sup_descubierta       0.00694 0.0000682    102.   0.         0.00680   0.00707

tidy(
  anova(modelo_log)
)

## # A tibble: 7 x 6
##   term                    df sumsq    meansq statistic p.value
##   <chr>                <int> <dbl>     <dbl>     <dbl>   <dbl>
## 1 log(rooms)               1 3724. 3724.       145610.       0
## 2 log(bathrooms)           1 1550. 1550.        60606.       0
## 3 log(surface_covered)     1 1381. 1381.        54013.       0
## 4 property_type            2  172.   86.2        3371.       0
## 5 barrios                  2 1481.  740.        28948.       0
## 6 sup_descubierta          1  265.  265.        10350.       0
## 7 Residuals            32123  821.    0.0256       NA       NA

De nuestro modelo confirmamos que todos los coeficientes resultan significativos.

Ahora comparamos nuestro nuevo modelo_log con nuestro modelo_sup_descubierta.

modelos = list(
  modelo_log = modelo_log,
  modelo_sup_descubierta = modelo_sup_descubierta
)

purrr::map_df(modelos, broom::glance, .id = "model")

## # A tibble: 2 x 13
##   model r.squared adj.r.squared   sigma statistic p.value    df  logLik     AIC
##   <chr>     <dbl>         <dbl>   <dbl>     <dbl>   <dbl> <dbl>   <dbl>   <dbl>
## 1 mode…     0.913         0.913 1.60e-1    41902.       0     8  1.33e4 -26605.
## 2 mode…     0.921         0.921 3.98e+4    41717.       0     9 -3.86e5 771799.
## # … with 4 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>,
## #   nobs <int>

Observamos que el mismo explica mejor la variabilidad \(Rˆ2\)

Ahora vamos a analizar el cumplimiento de los supuestos del modelo lineal

plot(modelo_log)

- Residuos vs valores predichos: parece existir una cierta estructura en los datos. La linealidad se respeta. Podemos observar un comportamiento heterocedástico aunque parece ser más uniforme que el del modelo_sup_descubierta. Se resaltan 3 observaciones (observaciones 3636,3709,7346). - Normal QQ plot: los extremos no se ajustan a la distribución teórica, tanto inferior izquierdo como superior derecho. Sin embargo parece que se comportan ligeramente mejor que el modelo_sup_descubierta (parecen estár mas suavizados los extremos). - Residual vs leverage:. se pueden apreciar 2 puntos con un leverage bastante alto y el grafico resalta 3 observaciones.

• Diagnóstico del modelo: El modelo creado no cumple con los supuestos del modelo lineal.Sin embargo y de observar los gráficos, parece que se comporta mejor que el modelo_sup_descubierta.

Selección del Modelo

Vamos a generar un nuevo modelo utilizando como variables \(log(rooms)\), \(log(bathrooms)\), \(log(surface_total)\), property_type y barrios y otro similar a nuestro modelo_log reemplazando la variable barrios por la variable l3.

Generación Nuevos Modelos

Generamos nuestro modelo modelo_log_surface_total que buscará inferir el resultado del \(log(price)\)

modelo_log_surface_total = lm(
  log(price) ~ log(rooms) + log(bathrooms) + log(surface_total) + property_type + barrios,
  data = train_ds
)
tidy_mlst <- tidy(modelo_log_surface_total, conf.int = TRUE)
tidy_mlst

## # A tibble: 8 x 7
##   term                 estimate std.error statistic   p.value conf.low conf.high
##   <chr>                   <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)           8.38      0.0139    603.    0.         8.35      8.40   
## 2 log(rooms)            0.00164   0.00275     0.594 5.53e-  1 -0.00376   0.00703
## 3 log(bathrooms)        0.0772    0.00297    26.0   2.52e-147  0.0714    0.0831 
## 4 log(surface_total)    0.935     0.00314   298.    0.         0.928     0.941  
## 5 property_typeDepart…  0.130     0.00566    23.0   9.67e-116  0.119     0.141  
## 6 property_typePH       0.0441    0.00599     7.36  1.90e- 13  0.0323    0.0558 
## 7 barriosprecio_bajo   -0.732     0.00258  -284.    0.        -0.737    -0.727  
## 8 barriosprecio_medio  -0.362     0.00207  -175.    0.        -0.366    -0.358

En nuestro modelo podemos observar como el coeficiente de \(log(rooms)\) parece no ser significativo (su p-valor = 0.5) y su intervalo de confianza incluye el 0. El resto de nuestros coeficientes parecen ser significativos.

Debido a que neustros coeficientes para la variable \(log(rooms)\) parece que no son significativos, vamos a analizar si la interacción de la variable si lo es para nuestro modelo.

tidy(
  anova(modelo_log_surface_total)
)

## # A tibble: 6 x 6
##   term                  df sumsq    meansq statistic p.value
##   <chr>              <int> <dbl>     <dbl>     <dbl>   <dbl>
## 1 log(rooms)             1 3724. 3724.       166879.       0
## 2 log(bathrooms)         1 1550. 1550.        69458.       0
## 3 log(surface_total)     1 1332. 1332.        59708.       0
## 4 property_type          2  263.  131.         5889.       0
## 5 barrios                2 1808.  904.        40523.       0
## 6 Residuals          32124  717.    0.0223       NA       NA

Y resulta que si lo es.

Analizamos el cumplimiento de los supuestos de nuestro modelo.

plot(modelo_log_surface_total)

Podemos ver como nuestro modelo se comporta de forma muy similar al modelo_log. Sin embargo hay una diferencia importante al analizar Residuos vs Leverage. Aca podemos observar como hay una especie de corte en el leverage entre los valores de 0.0007 y 0.0013 (aproximadamente).

Generamos nuestro modelo modelo_log_surface_total que buscará inferir el resultado del \(log(price)\)

modelo_log_l3 = lm(
  log(price) ~ log(rooms) + log(bathrooms) + log(surface_covered) + property_type + sup_descubierta + l3 ,
  data = train_ds
)
tidy_mll3 <- tidy(modelo_log_l3, conf.int = TRUE)
tidy_mll3

## # A tibble: 63 x 7
##    term                estimate std.error statistic   p.value conf.low conf.high
##    <chr>                  <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
##  1 (Intercept)          8.53    0.0265      322.    0.         8.48     8.58    
##  2 log(rooms)          -0.00924 0.00429      -2.15  3.12e-  2 -0.0177  -0.000834
##  3 log(bathrooms)       0.171   0.00431      39.7   0.         0.163    0.180   
##  4 log(surface_covere…  0.783   0.00506     155.    0.         0.773    0.793   
##  5 property_typeDepar…  0.198   0.00840      23.6   2.20e-122  0.182    0.215   
##  6 property_typePH      0.0476  0.00873       5.46  4.87e-  8  0.0305   0.0648  
##  7 sup_descubierta      0.00390 0.0000907    42.9   0.         0.00372  0.00407 
##  8 l3Agronomía          0.00514 0.0342        0.151 8.80e-  1 -0.0618   0.0721  
##  9 l3Almagro            0.00401 0.0166        0.242 8.09e-  1 -0.0286   0.0366  
## 10 l3Balvanera         -0.153   0.0176       -8.69  3.91e- 18 -0.188   -0.119   
## # … with 53 more rows

En este modelo, podemos observar el mismo problema que surgía al incorporar al analisis la variable l3. Sin embargo, el resto de nuestros coeficientes parecen ser estadisticamente significativos.

Debido a que neustros coeficientes para la variable l3 parece que no son significativos, vamos a analizar si la interacción de la variable si lo es para nuestro modelo.

tidy(
  anova(modelo_log_l3)
)

## # A tibble: 7 x 6
##   term                    df  sumsq    meansq statistic p.value
##   <chr>                <int>  <dbl>     <dbl>     <dbl>   <dbl>
## 1 log(rooms)               1 3724.  3724.        80828.       0
## 2 log(bathrooms)           1 1550.  1550.        33642.       0
## 3 log(surface_covered)     1 1381.  1381.        29982.       0
## 4 property_type            2  172.    86.2        1871.       0
## 5 sup_descubierta          1   81.3   81.3        1765.       0
## 6 l3                      56 1008.    18.0         391.       0
## 7 Residuals            32069 1477.     0.0461       NA       NA

Y resulta que si lo es.

Analizamos el cumplimiento de los supuestos de nuestro modelo.

plot(modelo_log_l3)

- Residuos vs valores predichos: parece existir una cierta estructura en los datos. La linealidad se respeta. Podemos observar un comportamiento heterocedástico (muy similar al modelo_log). Se resaltan 3 observaciones (observaciones 7346,7942, 28022). - Normal QQ plot: los extremos no se ajustan a la distribución teórica, tanto inferior izquierdo como superior derecho. Sin embargo parece haber cierto suavizado en el extremo superior derecho en comparación a otros modelos. - Residual vs leverage:. se pueden apreciar 3 puntos con un leverage bastante alto y el grafico resalta 3 observaciones (3636, 22246,26579).

• Diagnóstico del modelo: El modelo creado no cumple con los supuestos del modelo lineal. Parece que se comporta muy similar al modelo_log con una excepción al analizar el grafico de Residual vs leverage.

Evaluación de Modelos Seleccionados

Ahora vamos a elegir 2 de los modelos previamente desarrollados para evaluar y seleccionar, junto a nuestros 2 nuevos modelos, cual es el que mejor permite inferir el precio de una propiedad. Vamos a elegir los modelos que mejor expresan la variabilidad \(Rˆ2\). Estos son el modelo_log y modelo_sup_descubierta.

modelos = list(
  modelo_log = modelo_log,
  modelo_sup_descubierta = modelo_sup_descubierta,
  modelo_log_surface_total = modelo_log_surface_total,
  modelo_log_l3 = modelo_log_l3
)

purrr::map_df(modelos, broom::glance, .id = "model") %>%
  arrange(adj.r.squared)

## # A tibble: 4 x 13
##   model r.squared adj.r.squared   sigma statistic p.value    df  logLik     AIC
##   <chr>     <dbl>         <dbl>   <dbl>     <dbl>   <dbl> <dbl>   <dbl>   <dbl>
## 1 mode…     0.843         0.842 2.15e-1     2772.       0    62  3.88e3  -7638.
## 2 mode…     0.913         0.913 1.60e-1    41902.       0     8  1.33e4 -26605.
## 3 mode…     0.921         0.921 3.98e+4    41717.       0     9 -3.86e5 771799.
## 4 mode…     0.924         0.924 1.49e-1    55553.       0     7  1.55e4 -30987.
## # … with 4 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>,
## #   nobs <int>

Podemos observar que el modelo que mejor explica la variabilidad es el modelo modelo_log_surface_total. Observamos el valor del \(Rˆ2\) ajustado.

Predicción

Cargamos el dataset de testing y le aplicamos nuestras transformaciones.

propiedades_test <- read_csv(here("/ds/ar_properties_test.csv"))

## Parsed with column specification:
## cols(
##   id = col_character(),
##   l3 = col_character(),
##   rooms = col_double(),
##   bathrooms = col_double(),
##   surface_total = col_double(),
##   surface_covered = col_double(),
##   price = col_double(),
##   property_type = col_character()
## )

propiedades_test = propiedades_test %>%
  mutate(
    precio_m2 = price / surface_total,
    barrios = factor(
      case_when(
        precio_m2 <= quantile(precio_m2)[2] ~ "precio_bajo",
        precio_m2 > quantile(precio_m2)[2] & precio_m2 <= quantile(precio_m2)[4] ~ "precio_medio",
        precio_m2 > quantile(precio_m2)[4] ~ "precio_alto"
      )
    ),
    sup_descubierta = surface_total - surface_covered
  )

propiedades_test

## # A tibble: 13,772 x 11
##    id    l3    rooms bathrooms surface_total surface_covered  price
##    <chr> <chr> <dbl>     <dbl>         <dbl>           <dbl>  <dbl>
##  1 Afdc… Vele…     3         2            95              69 199900
##  2 ESzy… Nuñez     1         1            44              38 147000
##  3 R7Is… Alma…     1         1            40              37  77000
##  4 MFng… Alma…     1         1            40              37  92943
##  5 j3TP… Pale…     1         1            32              30 125000
##  6 jtR+… Pale…     1         1            45              35 150000
##  7 UOCj… Pale…     1         1            33              32 110000
##  8 3mG1… Parq…     1         1            70              70 225000
##  9 Nmld… Alma…     2         1            48              48 129900
## 10 LmF0… Reti…     2         1            58              58 158000
## # … with 13,762 more rows, and 4 more variables: property_type <chr>,
## #   precio_m2 <dbl>, barrios <fct>, sup_descubierta <dbl>

Evaluamos en Training

Ahora vamos a evaluar nuestros modelos en training y evaluar sus predicciones en testing. Para esto es necesario hacer un análisis distinto para los modelos que utilizan \(log\) que para el modelo_sup_descubierta.

modelos_log = list(
  modelo_log = modelo_log,
  modelo_log_l3 = modelo_log_l3,
  modelo_log_surface_total = modelo_log_surface_total
)

predicciones_training_log = map(.x = modelos_log, .f = augment)

map_dfr(
  .x = predicciones_training_log, 
  .f = rmse, 
  truth = exp(`log(price)`), 
  estimate = exp(.fitted), 
  .id="modelo"
) %>% arrange(.estimate)

## # A tibble: 3 x 4
##   modelo                   .metric .estimator .estimate
##   <chr>                    <chr>   <chr>          <dbl>
## 1 modelo_log_surface_total rmse    standard      44762.
## 2 modelo_log               rmse    standard      47444.
## 3 modelo_log_l3            rmse    standard      62055.

eval_train_sd = augment(modelo_sup_descubierta)

rmse(
  data = eval_train_sd,
  truth = price,
  estimate = .fitted
)

## # A tibble: 1 x 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard      39757.

Analizando los RMSE de nuestros modelos en training, podemos observar que el modelo con el menor valor es nuestro modelo_sup_descubierta seguido (en orden) por los modelos modelo_log_surface_total, modelo_log, modelo_log_l3.

Evaluamos en Testing

Ahora vamos a evaluar nuestros modelos en testing.

# Aplicamos la función augment a los 4 modelos con el set de testing
predicciones_testing = map(
  .x = modelos_log, 
  .f = augment,
  newdata = propiedades_test
) 

# Obtenemos el RMSE para los 4 modelos
map_dfr(
  .x = predicciones_testing, 
  .f = rmse, 
  truth = price, 
  estimate = exp(.fitted), 
  .id="modelo"
) %>% 
  arrange(.estimate)

## # A tibble: 3 x 4
##   modelo                   .metric .estimator .estimate
##   <chr>                    <chr>   <chr>          <dbl>
## 1 modelo_log_surface_total rmse    standard      43832.
## 2 modelo_log               rmse    standard      46153.
## 3 modelo_log_l3            rmse    standard      60873.

pred_log = augment(
  modelo_sup_descubierta, 
  newdata=propiedades_test
) 
pred_log

## # A tibble: 13,772 x 13
##    id    l3    rooms bathrooms surface_total surface_covered  price
##    <chr> <chr> <dbl>     <dbl>         <dbl>           <dbl>  <dbl>
##  1 Afdc… Vele…     3         2            95              69 199900
##  2 ESzy… Nuñez     1         1            44              38 147000
##  3 R7Is… Alma…     1         1            40              37  77000
##  4 MFng… Alma…     1         1            40              37  92943
##  5 j3TP… Pale…     1         1            32              30 125000
##  6 jtR+… Pale…     1         1            45              35 150000
##  7 UOCj… Pale…     1         1            33              32 110000
##  8 3mG1… Parq…     1         1            70              70 225000
##  9 Nmld… Alma…     2         1            48              48 129900
## 10 LmF0… Reti…     2         1            58              58 158000
## # … with 13,762 more rows, and 6 more variables: property_type <chr>,
## #   precio_m2 <dbl>, barrios <fct>, sup_descubierta <dbl>, .fitted <dbl>,
## #   .resid <dbl>

rmse(
  data = pred_log, 
  truth = price, 
  estimate = .fitted
)

## # A tibble: 1 x 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard      38975.

Volvemos a observar como el modelo con el menor valor de RMSE es nuestro modelo_sup_descubierta.

Conclusion

En lineas generales no podemos afirmar que nuestros modelos cumplen los supuestos del modelo lineal.

Sin embargo, podemos utilizar nuestros modelos para predecir los valores del precio de las propiedades. Ahora bien, para elegir el mejor de los modelos propuestos, vamos a fundamentar nuestra decisión en base a la evaluación de las 2 métricas que estamos observando: - \(Rˆ2\) - RMSE

Los modelos que mejor explican la variabilidad medida por \(Rˆ2\) ajustado son el modelo_log, modelo_sup_descubierta y modelo_log_surface_total. Todos estos modelos se encuentran con valores superiores a 0.91.

Sin embargo, nosotros queremos predecir nuevos datos y para ello es importante medir el RMSE para evaluar el error en la predicción. Al analizar esta metrica, observamos un salto importante entre nuestro modelo_sup_descubierta (38974.87) y sus seguidores modelo_log_surface_total (43832.07) y modelo_log (46152.97).

Por ende seleccionaría como mi mejor modelo el modelo modelo_sup_descubierta.